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A  THREE-DIMENSIONAL  CALCULATION  OF 
THE  VENTING  PROCESS  FOR  AN 
ON-AXIS  EXPLOSION  IN  THE  MX  TRENCH 


ABSTRACT 

Using  the  TRIDORF  code,  we  performed  a  three-dimensional  calculation  to 
study  the  venting  of  energy  through  the  roof  of  the  proposed  MX  missile  emplacement 
trench  for  an  on-axis  detonation  of  1  Mt.  A  two-dimensional  calculation  to  27  ns  by  Science 
Applications,  Incorporated,  gave  us  the  initial  conditions;  our  calculation  extended  to 
238  /iS,  through  the  period  in  which  venting  and  asymmetrical  tunnel  expansion  are  the  ma¬ 
jor  energy-transfer  processes  in  the  trench.  These  two  calculations  provide  the  best  available 
starting  conditions  for  later  1-1/2-  and  2-dimensional  calculations  of  the  downstream  flow. 


INTRODUCTION 


One  of  the  options  considered  by  the  Air  Force 
for  land-based  weapon  systems  is  the  MX  trench. 
To  assess  the  degree  of  protection  afforded  the  mis¬ 
sile  by  blast  doors  and  ribs  in  the  trench  walls, 
calculations  have  been  made  with  quasi-one- 
dimensional  and  two-dimensional  codes  of  the 
trench  environment  during  an  enemy  attack  for 
near-misses,  direct  hits,  and  in-trench  explosions.*'^ 
Because  the  problem  is  inherently  three- 
dimensional,  various  leakage  models  have  been 
used  to  simulate  venting  through  the  tunnel  roof.  In 
terms  of  machine  size  and  running  time,  quasi-one- 
dimensional  codes  are  the  most  practical  means  of 
carrying  out  calculations  over  long  distances  and  to 
very  late  times.  For  shorter  distances  and  shorter 
times,  to  approximate  the  three-dimensional  con¬ 
figuration,  two-dimensional  calculations  have  been 
made  along  the  plane  of  symmetry  of  the  problem 
and  for  a  cross-sectional  plane  through  the  source. 
However,  from  the  onset  of  venting  through  the 


period  in  which  venting  is  the  major  energy-transfer 
process,  when  the  phenomenology  of  interest  is 
truly  three-dimensional,  three-dimensional  calcula¬ 
tions  are  needed. 

To  investigate  trench  conditions  during  the 
restricted  period  in  which  roof  venting  is  expected 
to  be  the  dominant  transfer  mechanism,  we  per¬ 
formed  a  three-dimensional  calculation  using  the 
TRIDORF  version  of  W.  E.  Johnson’s  TRIOIL 
code.  TRIDORF,^  a  two-material,  three- 
dimensional,  Eulerian  hydrodynamic  code,  has  a 
rigid-plastic  strength  model  and  radiation  diffusion 
capability.  The  calculation  run  as  a  purely 
hydrodynamic  problem  extended  from  about  27  //s, 
when  the  problem  is  just  becoming  three- 
dimensional,  to  about  260  ^s.  This  report  describes 
the  code  and  calculational  setup,  gives  results  of  our 
calculations  to  238  ^s,  and  compares  our  results  at 
260  ^s  with  those  of  quasi-one-dimensional  and 
two-dimensional  calculations. 


CODE  AND  CALCULATIONAL  SETUP 


The  original  problem  involved  a  1-Mt  source 
detonated  on  the  axis  of  a  tunnel  4.2  m  in  diameter. 
The  tunnel  walls  were  assumed  to  be  iron-fiber- 
filled  concrete  0.275  m  thick.  The  surrounding 
medium  was  alluvium  and  the  overburden  thickness 
1.495  m. 


We  generated  the  three-dimensional  problem 
by  linking  to  a  two-dimensional  calculation  per¬ 
formed  by  SAl  (Science  Applications,  Incor¬ 
porated)  using  the  DRAGNET  and  RADOIL 
codes.  The  link  was  made  at  27.06  fis,  when  the 
problem  was  still  two-dimensional.  The  TRIDORF 


I 


generator  spun  the  two-dimensional  problem  about 
its  axis  to  form  half  the  cylindrical  tunnel  and  also 
added  air  above  the  tunnel  and  more  alluvium 
below,  forming  the  three-dimensional  problem  il¬ 
lustrated  in  Fig.  1. 

Figures  2  through  4  show  the  zoning  for  the 
TRIDORF  problem  along  with  the  corresponding 
material  maps  at  link  time.  In  all,  the  calculation 
used  36,000  zones;  30  zones  (JMAX)  covering 
87.37  m  from  the  source  along  the  tunnel  length,  50 
zones  (IMAX)  extending  41.18  m  in  the  direction 
perpendicular  to  both  the  tunnel  axis  and  the 
ground  surface,  and  24  zones  (KMAX)  extending 
8.44  m  in  the  third  direction.  The  two-dimensional 
RADOIL  problem  to  which  TRIDORF  was  linked 
had  61  zones  along  the  tunnel  axis  and  137  in  the 
radial  direction.  Because  of  storage  limitations  in 
the  CDC  7600  computer,  the  full  length  of  the  tun¬ 
nel  in  the  RADOIL  problem  was  not  included  in  the 
TRIDORF  calculation  but  only  the  87  m  (55  of  61 
zones)  needed  to  carry  it  past  the  time  when  venting 
is  expected  to  dominate.  Large  “sink”  zones  were 
placed  at  the  edges  of  the  problem  except  at  the 
plane  of  symmetry  and  the  tunnel  ends.  To  allow  a 
better  resolved  calculation  of  the  venting,  more 
zones  were  used  on  the  external-air  side  of  the 


problem  above  the  trench  than  below  it.  Transmit- 
tive  boundary  conditions  allowing  the  passage  of 
momentum  and  energy  over  the  grid  boundaries 
were  applied  except  across  the  plane  of  symmetry. 

To  maintain  consistency  with  the  two- 
dimensional  SAI  calculation,  the  equations  of  state 
in  TRIDORF  were  made  identical  with  those  in 
their  version  of  RADOIL.  Because  both  are  intrin¬ 
sically  two-material  codes,  the  tunnel  wall  material 
and  tunnel  air  were  labeled  as  one  material  (dot) 
distinguished  only  by  density.  When  the  density  of 
the  dot  material  was  less  than  0.005  g/cm^,  it  was 
treated  as  air  and  described  as  an  ideal  gas  with  a  7 
of  1.5.  For  greater  densities  it  was  treated  as  con¬ 
crete  using  a  modiOed  form  of  the  Tillotson  equa¬ 
tion.  The  constants  for  the  concrete  are  given  in 
Table  1  along  with  those  for  the  alluvium  (X 
material)  outside  the  trench  wall. 

The  modified  Tillotson  equation  of  state  for 
the  dot  and  X  materials  was  as  follows.  Let  pp  be  the 
normal  density,  p  the  density  of  the  material  in  the 
zone,  P  the  pressure,  and  I  the  specific  internal 
energy.  The  constant  A  is  the  bulk  modulus;  G  is 
the  zero-pressure  Griineisen  constant.  Iq,  b,  B,  «, 
and  li  are  fitting  constants.  The  quantity  1^  is  the 
specific  internal  energy  of  the  material  at  the 


FIG.  1.  MX  configuration  at  three-dimensional  link  time  (27  ps). 
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FIG.  2.  Symmetry  plane:  (a)  X-Y  zoning;  (b)  material  map. 
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FIG.  3.  Cross-sectional  plane:  (a)  Z-X  zoning;  (b)  material  map. 
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FIG.  4.  Longitudinal  plane:  (a)  Y-Z  zoning;  (b)  material  map  for  40-m  area  closest  to  source. 


vaporization  temperature;  ij  is  equal  to  plus  the 
energy  needed  to  vaporize  the  material  from  the  liq¬ 
uid  state.  Thus; 


n  =  P/Pq 
P  =  p/pQ  -  1 


ForT?>l: 

P(p,  I)  =  Gpl  + 

For  7?  <  1  and  1  <  1^; 

P>,I)  =  GpI  +  Ap 

For7?<  1  and  I  >1': 

P.  (p,  I)  =  alp  +  (  f-  .  - 

UK'o’'  )J*' 

t  '“'o'" -"V 


For77<landIj<I<i;: 

PC.,,-  p.cp.,, 


TABLE  1.  TOlotaon  constants  for  concrete  and  alluvium. 


Conttmt 

Concrete 

Allimum 

P0.«/cin^ 

2.5 

1.76 

M 

0.5 

0.5 

Eo.kJ/g 

ILO 

0.00201 

b 

IJ 

IJ 

A,kJ/cin^ 

4j6 

23A 

E^,U/g 

1.91 

3.50 

E^.U/g 

3.20 

182) 

a 

52) 

52) 

9 

5j0 

5.0 

B,kJ/g 

23.5 

10.2 

For  concrete  vapor  of  density  less  than  5  X  I0~\  a 
gamma-law  equation  of  state  for  air  was  used  with 
(7  -  I)  =  0.5. 

The  total  energy  in  the  problem  at  link  time  (27 
ps)  was  1.025  X  10'^  J  (0.25  Mt),  representing  the 
quarter  of  the  total  1-Mt  yield  that  was  included  in 
the  grid  (see  geometry  in  Fig.  1).  At  175  ps  after 
detonation,  1%  of  the  total  energy  had  flowed  over 
the  grid  boundary.  At  250  ps,  the  end  of  the  calcula¬ 
tion,  the  loss  had  increased  to  5.8%  of  the  total 
energy.  The  maximum  amount  of  radiation  energy 
in  any  zone  at  link  time  was  about  2%  of  the  total 
internal  energy  in  that  zone;  hence,  the  radiation 
option  in  the  TRIDORF  code  was  not  used.  The 
calculation  from  27  to  250  ps  required  some  2  h  of 
CDC  7600  computer  time. 


CALCULATIONAL  RESULTS  TO  238  ms 


Supporting  the  following  description  are  plots 
of  variables  made  with  the  TENPLT  code^  at  and 
after  link  time.  The  calculational  data  allow  detailed 
analysis  of  conditions  in  the  trench  during  the  time 
in  which  venting  and  asymmetrical  tunnel  expan¬ 
sion  are  the  major  energy-transfer  mechanisms. 

Figure  5  is  a  vector  plot  of  velocities  at  link 
time.  The  apparent  rotation  of  the  vectors  back 
toward  the  tunnel  axis  in  the  region  Y  =40  m  occurs 
because  the  vectors  are  plotted  at  locations  along 
straight  lines  parallel  to  the  problem  axis  rather 
than  along  isovelocity  lines.  Figure  6  is  an  isometric 
plot  of  the  tunnel  pressure  at  link  time.  The  pressure 


front  traveling  radially  has  not  yet  broken  through 
the  alluvium  overburden  into  the  air.  The  accom¬ 
panying  geometrical  view  serves  merely  to  orient  the 
reader,  and  no  distinction  is  made  between  the  wall 
material  and  the  alluvium.  The  shaded  portion  in¬ 
cludes  both. 

Figure  7  shows  plots  of  the  pressure  proiile 
along  the  tunnel  axis  at  several  times  during  the 
calculation.  Figure  8  is  a  semilog  plot  of  these 
pressure  profiles  for  the  range  beyond  30  m.  Here 
one  can  see  at  about  60  m  the  diffusion  front  es¬ 
tablished  at  early  times  in  the  two-dimensional 
calculation,  which  was  able  to  move  no  farther 


down  the  tunnel  because  the  coupling  of  most  of  its 
energy  into  the  wall  left  little  radiant  energy  to  drive 
it.  By  about  80  ns  this  front  has  collapsed  and  has 
been  overtaken  by  hydrodynamic  pressure  moving 
down  the  tunnel.  Figure  9  shows  the  density  profile 
along  the  tunnel  axis  at  28.7  /is  and  at  238  /is.  The 
location  of  the  diffusion  front  is  again  apparent  in 
the  density  plot  at  the  earlier  time. 

Figure  10  shows  the  magnitudes  of  the  Y- 
directed  velocities  along  the  tunnel  at  the  axis. 
Figure  1 1  is  the  material  map  at  238  /is.  The  venting 
process  is  very  apparent  here.  Figure  12  is  a  blowup 
of  the  material  plot  in  the  region  of  the  venting 
showing  velocity  vectors  in  that  region.  The  vector 
arrows  are  not  intended  to  show  velocity 
magnitudes  but  indicate  the  direction  and  location 
of  the  stagnation  region  at  the  trench  bottom. 
(Single-headed  arrows  are  velocities  larger  than  the 
scale  maximum.)  Similar  material  and  velocity  vec¬ 
tors  are  shown  in  Fig.  13  for  the  tunnel  region 
beyond  40  m  from  the  source. 

Figures  14  through  17  show  isometric  plots  of 
density  and  pressure,  along  with  corresponding 
geometrical  views.  By  238  /is  the  trench  top  has  been 
unzipped  to  about  30  m.  This  is  also  the  approx¬ 
imate  location  of  the  pressure  peak  along  the  tunnel 
axis  at  this  time,  as  can  be  seen  in  Fig.  18.  Figures  19 
through  22  show  the  pressure  profiles  along  vertical 
lines  through  the  tunnel  axis  from  the  trench  bot¬ 


tom  up  through  the  overburden  at  increasing  dis¬ 
tances  down  the  tunnel.  At  238  ns  the  problem  is 
still  two-dimensional  beyond  50  m.  Figures  23 
through  34  are  cross-sectional  material  maps,  vector 
velocity  plots,  and  isometric  plots  of  density  and 
pressure  for  distances  of  0.7,  28,  and  50  m  from  the 
source.  Figure  35  shows  the  trench  configuration  on 
the  symmetry  plane  of  the  problem  at  238  /is.  Tun¬ 
nel  explosion  has  extended  out  to  about  46  m  from 
the  source,  but  there  is  still  material  of  significant 
density  (greater  than  1  g/cm^)  in  the  overburden 
back  to  about  29  m.  The  wall  vapor-air  interface  as 
defined  by  the  boundary  density  of  0.005  g/cm^  has 
moved  down  the  tunnel  to  about  60  m. 

We  calculated  the  fraction  of  total  energy 
remaining  in  the  original  tunnel  volume  using  the 
post-processer  code  SAMPLER. **  Values  as  a  func¬ 
tion  of  time  are  shown  in  Fig.  36.  The  most  rapid 
loss  at  early  times  results  from  radiant  coupling  of 
the  energy  to  the  wall,  leading  to  rapid  venting  while 
the  radiation  front  remains  stationary  in  the  tunnel. 
Venting  slows  as  the  hydrodynamic  pressure  front 
moves  down  the  tunnel  and  overtakes  the  stagnant 
radiation  front,  producing  the  knee  in  the  curve  at 
about  60-70  /is.  Venting  due  to  the  hydrodynamic 
pressure  continues  at  a  slower  rate;  at  260 /is,  about 
9%  of  the  total  energy  remains  in  the  original  trench 
volume. 


COMPARISON  WITH  ONE-  AND 
TWO-DIMENSIONAL  CALCULATIONS 


SAI  had  performed  a  two-dimensional  calcula¬ 
tion  in  the  symmetry  plane  of  the  three-dimensional 
TRIDORF  calculation  and  also  in  the  cross- 
sectional  plane  through  the  source.^  Figure  37  com¬ 
pares  peak  pressures  as  a  function  of  time  for  these 
calculations.  At  40  /is  stagnation  occurs  on-axis, 
producing  a  peak  in  the  profile.  In  a  two- 
dimensional  code  this  pressure  should  be  better 
calculated  in  transverse  than  in  longitudinal 
geometry  and,  indeed,  we  find  excellent  agreement 
between  the  transverse  two-dimensional  and 
TRIDORF  calculations.  In  general,  the  peak 
pressures  by  TRIDORF  are  30-35%  lower  than  by 
the  two-dimensional  calculations.  This  is  to  be  ex¬ 
pected  because  the  three-dimensional  calculation 
has  the  additional  divergence  term.  The  two- 
dimensional  calculations  are  thus  somewhat  conser¬ 


vative.  Figures  38  through  40,  showing  density, 
pressure,  and  velocity  profiles  at  260  /is,  also  exhibit 
this  tendency. 

Figure  41  shows  pressure  vs  distance  down  the 
tunnel.  Comparison  is  made  between  the  SAI  quasi- 
one-dimensional  (STEALTH)  calculation  and  the 
on-axis  values  for  the  TRIDORF  calculation  at  250 
/IS.  The  peak  pressure  has  advanced  about  100  m 
farther  downstream  in  the  one-dimensional  calcula¬ 
tion  than  in  the  three-dimensional  calculation,  but 
the  values  of  peak  pressure  for  the  two  calculations 
are  almost  identical  at  250  /is.  Different  venting 
models  might  be  tried  for  the  one-dimensional 
calculation  to  reduce  the  discrepancy.  This  would 
be  worthwhile  because  the  one-dimensional  calcula¬ 
tion  is  the  only  practical  one  to  pursue  at  very  late 
times. 
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FIG.  5.  Velocities  on  symmetry  pluie  at  link  time. 
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FIG.  6.  Pressure  isometric  plot  at  link  time. 
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FIG.  7.  Pressure  along  tunnel  axis. 
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FIG.  9.  Density  profile  along  tunnel  axis. 
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FIG.  12.  Expanded  material  plot  with  velocity  vectors  on  symmetry  plane  at  238  ms  (source  end). 
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FIG.  17.  Isometric  plot  of  pressure  on  symmetry  piane  at  238  fts  (looking  from  downstream  end). 


18 


Pressure  —  Mb 


FIG.  18.  Pressure  profile  along  tunnel  axis  at  238  ns. 
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FIG.  19.  Pressure  profile  across  tunnel  10  m  from  source. 
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FIG.  21.  Pressure  profile  across  tunnel  25  m  from  source. 
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FIG.  22.  Pressure  profile  across  tunnel  SO  m  from  source. 
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Time:  238  fis 


FIG.  23.  Material  map  with  vdocity  vectors  0.7  m  from  source. 
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Kploded  velocity  vectors  0.7  m  from  source. 
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FIG.  25.  Isometric  density  plot  0.7  m  from  source. 
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FIG.  26.  Isometric  pressure  plot  0.7  m  from  source. 
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FIG.  27.  Material  map  with  velocity  vectors  28  m  from  source. 
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FIG.  28.  Isometric  density  plot  28  m  from  source. 
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FIG.  29.  Isometric  density  plot  28  m  from  source. 
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FIG.  30.  Isometric  pressure  plot  28  m  from  source. 
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FIG.  31.  Material  map  with  velocity  vectors  SO  m  from  source. 
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FIG.  32.  Isometric  density  plot  50  m  from  source. 
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FIG.  33.  Isometric  density  plot  50  m  from  source. 
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FIG.  34.  Isometric  pressure  plot  50  m  from  source. 
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FIG.  35.  Geometric  'xmfiguratioa  oo  symmetry  plane  at  238  ms. 


36 


ooixra/3) 


Time  —  ns 

FIG.  36.  Fraction  of  total  energy  in  original  tunnel  volnme  as  a  function  of  time.  Data  points  are  taken  at  edit 
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FIG.  37.  Peak  pressure  on  tunnel  axis  as  a  function  of  time  comparing  results  of  three-dimensional  calculations 
with  those  of  transverse  and  longitudinal  two-dimensionai  calculations. 
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FIG.  38.  Comparison  of  three-dimensional  and  longitudinal  two-dimensional  calculations:  density. 
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FIG.  39.  Comparison  of  three-dimensional  and  longitudinal  two-dimensional  calculations:  pressure. 
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FIG.  40.  Comparison  of  three-dimensional  and  longitudinal  two-dimensional  calculations:  velocity  along  tunnel 
axis. 


FIG.  41.  Comparison  of  presMire  profiles  by 
three-<fimeasional  and  one-dfaneosional  calcula- 
dons  at  250  fi%. 
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CONCLUSIONS 


The  three-dimensional  TRIDORF  code  has 
been  used  si  ■  ?ssfully  to  calculate  the  venting  that 
would  re  in  a  l-Mt  detonation  on-axis  in  the 
MX  treni.  ;is  study  demonstrates  both  the  need 
for  and  the  practicality  of  using  a  three-dimensional 
code  for  calculations  in  which  the  phenomenology 
of  interest  is  truly  three-dimensional.  Traditionally 
it  has  been  argued  that  one  should  perform  several 
crude  one-dimensional  calculations  followed  by  one 
or  two  selected  two-dimensional  calculations  and 
only  then,  if  necessary,  by  a  three-dimensional 
calculation.  Here,  however,  the  phenomenology 


allows  an  early-time  two-dimensional  study  but  re¬ 
quires  a  three-dimensional  calculation  to  span  the 
period  of  venting.  At  the  end  of  this  time,  these  two 
calculations  can  provide  satisfactory  initial  condi¬ 
tions  for  later  1-1/2-  and  2-dimensional  studies  of 
the  downstream  flow. 

The  three-dimensional  calculation  was  run  in  a 
straightforward  manner  without  the  need  for 
postulating  special  models  to  account  for  the 
venting  process.  The  36,000-zone  problem  required 
just  2  h  of  CDC  7600  computer  time. 


REFERENCES 

1.  A.  L.  Kuhl,  Estimates  of  Airblast  Environments  Inside  the  MX  Buried  Trench,  TRW  Defense  and  Space 
Systems  Group,  Redondo  Beach,  CA,  Report  76.435.9-72  (1976). 

2.  C.  E.  Needham  and  C.  Westmoreland,  Airblast  Calculations  for  Advanced  Missile  Systems  ( MX  Support), 
Air  Force  Weapons  Laboratory,  Albuquerque,  NM,  Report  AFWL-TR-75-297  (1976).  (Distribution 
limited  to  U.S.  Government  agencies  only;  AD-B-009  6502,  Defense  Documentation  Center,  Alexandria, 
VA.) 

3.  R.  N.  Schlaug  and  B.  F.  Chambers  III,  Theoretical  Calculations  in  Support  of  In-Trench  Environment 
Definition  for  the  MX-Concept  Validation  Phase,  Defense  Nuclear  Agency,  Report  DNA-4655F-2,  Vol.  2 
(1978). 

4.  W.  E.  Johnson,  TRIDORF:  A  Two-Material  3-D  Cartesian  Hydrodynamics  Code  Coupled  with  a  Rigid 
Plastic  Strength  Model,  Computer  Code  Consultants,  Incorporated,  Los  Alamos,  NM,  Report  CCC-976 
(1976). 

5.  D.  E.  Burton  and  C.  M.  Snell,  User’s  Guide  to  TENPLT  (TENSOR  Graphics  Code)  and  EOSLIB  AND 
XMUPGENS  (Auxiliary  Data  Codes),  U.  S.  Army  Engineer  Waterways  Experimental  Station  Explosives 
Evacuation  Research  Laboratory,  Livermore,  CA,  Miscellaneous  Paper  E-74-1  (1974). 


eSM/ckg 


U.S.  Government  Printing  Office:  1980/2  -  689-102/8753 


41 


